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Abstract.  A  shallow-water  model  was  used  to  understand 
model  error  induced  by  non-Gaussian  wind  uncertainty.  Al¬ 
though  the  model  was  simple,  it  described  a  generic  sys¬ 
tem  with  many  degrees  of  freedom  randomized  by  external 
noise.  The  study  focused  on  the  nontrivial  collective  behav¬ 
ior  of  finite-amplitude  perturbations  on  different  scales  and 
their  influence  on  model  predictability.  The  error  growth 
strongly  depended  on  the  intensity  and  degree  of  spatial  in¬ 
homogeneity  of  wind  perturbations.  For  moderate  but  highly 
inhomogeneous  winds,  the  error  grew  as  a  power  law.  This 
behavior  was  a  consequence  of  varying  local  characteris¬ 
tic  exponents  and  nonlinear  interactions  between  different 
scales.  Coherent  growth  of  perturbations  was  obtained  for 
different  scales  at  various  stages  of  error  evolution.  For  the 
nonlinear  stage,  statistics  of  prediction  error  could  be  ap¬ 
proximated  by  a  Weibull  distribution.  An  approach  based 
on  the  Kullback-Leibler  distance  (the  relative  entropy)  and 
probability-weighted  moments  was  developed  for  identifi¬ 
cation  of  Weibull  statistics.  Bifurcations  of  the  variance, 
skewness  and  kurtosis  of  the  irreversible  predictability  time 
(a  measure  of  model  prediction  skill)  were  detected  when 
the  accepted  prediction  accuracy  (tolerance)  exceeded  some 
threshold. 


1  Introduction 

When  circulation  is  simulated  by  a  fine  resolution  regional 
model  in  an  area  with  open  boundaries,  the  circulation  dy¬ 
namics  often  depends  crucially  on  specified  open  bound¬ 
ary  conditions,  wind  forcing  and  sub-scale  parameteriza- 
tions.  For  atmospheric  predictability,  it  is  generally  assumed 
that  the  model  forecasting  is  most  sensitive  to  uncertainty 
of  initial  conditions.  However,  for  oceanic  predictability  in 
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marginal  seas  with  lengthy  coastal  zones,  impact  of  uncer¬ 
tainty  of  external  forcing  and/or  subgrid  parameterizations, 
may  be  as  significant  as  errors  in  initial  conditions.  The 
model  blow-up  often  sets  in  much  faster  than  the  model  loses 
its  predictability  due  to  errors  in  initial  conditions  (Jiang 
and  Malanotte-Rizzoli,  1999;  Boffetta  et  ah,  2000;  Bogden, 
2001;  Auclair  et  ah,  2003,  among  others). 

Starting  from  the  pioneering  work  of  Lorenz  (1963),  the 
dynamics  of  the  prediction  error  (PE)  due  to  uncertainty  in 
initial  conditions  has  been  deeply  investigated  in  the  theoret¬ 
ical  and  numerical  studies.  The  dynamics  of  model-related 
errors  has  been  paid  much  less  attention  to,  probably  due 
to  the  large  variety  of  possible  modeling  errors.  Although 
some  general  trends  have  emerged  (Capotondi  and  Holland, 
1997;  Chu  et  ah,  1999;  Orrell  et  ah,  2001;  Vannitsem  and 
Toth,  2002;  Nicolis,  2003,  2004,  among  others),  more  refined 
theoretical  investigations  and  additional  experiments  with  an 
hierarchy  of  ocean  models  of  different  levels  of  complexity 
are  necessary  to  get  a  more  general  view  of  the  impact  of 
model-related  error  (particularly  finite-amplitude)  on  model 
predictability. 

In  the  present  paper,  the  dynamics  of  a  model-related  error 
(hereafter,  prediction  error  (PE))  generated  by  uncertainty  in 
wind  forcing  and  its  impact  on  predictability  are  studied  in 
the  context  of  stochastic  model  stability  (stability  answered 
in  terms  of  probabilistic  measures,  such  as  expected  values 
or  distribution  functions  (Freidlin  and  Wentzell,  1998))  for  a 
simplified  regional  model  destabilized  near  an  unstable  equi¬ 
librium  state  (an  unstable  fixed  point  in  model  phase  space) 
by  stochastic  wind. 

In  general,  the  stochastic  stability  and  predictability  dif¬ 
fer  from  one  another.  However,  if  a  time  scale  quantifies  the 
model  predictability,  and  if  this  scale  indicates  the  time  when 
the  forecast  uncertainty  exceeds  some  boundary  or  when  in¬ 
formation  on  the  initial  condition  is  lost,  the  stochastic  sta¬ 
bility  and  predictability  are  interchangeable.  Since  such  time 
scales  are  widely  used  in  meteorology  (see,  for  example. 
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Toth,  1991)  and  oceanography  (Robinson  et  al.,  1996),  the 
stochastic  stability  concept  seems  to  be  a  useful  tool  for  the 
predictability  analysis  of  large  ocean  models,  and  “the  irre¬ 
versible  predictability  time”  (a  prediction  time  scale)  (Chu  et 
al.,  2002)  is  a  quantitative  measure  of  model  predictability. 

Numerical  computations  discussed  below  can  have  a 
twofold  interpretation.  First,  the  obtained  results  can  be  in¬ 
terpreted  as  stability  of  a  circulation  regime  relative  to  a  high- 
frequency  component  of  wind  forcing.  Second,  if  we  pa¬ 
rameterize  uncertainty  of  wind  forcing  as  a  stochastic  noise, 
the  computation  results  are  interpretable  in  the  predictabil¬ 
ity  context.  Therefore,  hereafter,  perturbations  excited  by 
stochastic  wind  in  an  ocean  basin  will  also  be  called  a  pre¬ 
diction  error. 

The  principal  motivation  of  the  proposed  study  is  formu¬ 
lated  as  follows. 

First,  it  is  well-known  from  numerical  modeling  that  wind 
is  one  of  the  main  energy  sources  of  ocean  currents  and  that, 
naturally,  prediction  errors  associated  with  wind  forcing  un¬ 
certainty  may  grow  quickly  in  numerical  models  (Robinson 
et  al.,  1996;  Berloff  and  McWilliams,  1999;  Chu  et  al.,  1999; 
Sura  et  al.,  2001;  Burillo  et  al.,  2002,  and  others).  Behav¬ 
ior  of  such  model-related  errors  are  highly  dependent  on  the 
properties  of  the  underlying  dynamical  regime  [attractor]  re¬ 
produced  by  the  numerical  model  and  statistics  of  wind  un¬ 
certainty.  Circulation  patterns  with  oscillations  near  quasi¬ 
equilibrium  states  and  transition  dynamics  between  them 
are  typical  for  many  marginal  seas,  such  as  the  Black  Sea 
(Stanev  and  Staneva,  2000),  and  for  large-scale  jet-like  cur¬ 
rents  like  Kuroshio  (Masuda  et  al.,  1999).  The  proposed 
study  focuses  on  the  case  of  evolution  of  finite-amplitude 
non-Gaussian  perturbations  induced  by  stochastic  wind  error 
when  the  ocean  is  near  a  quasi-stable  equilibrium  state  (even 
small  perturbations  can  destroy  such  a  state  and  stimulate  a 
transition  to  another  one).  Therefore,  the  obtained  results  are 
important  for  understanding  the  regional  model  predictabil¬ 
ity  (Robinson  et  al.,  1996)  when  local  attractor  features  deter¬ 
mine  the  phase-spatial  organization  of  the  local  error  growth 
rate. 

Second,  in  ocean  models,  unresolved  dynamics  is  often 
represented  in  terms  of  random  forcing.  For  example,  impact 
of  mesoscale  eddies  on  large-scale  currents  can  be  approx¬ 
imated  as  a  space-time  correlated,  random-forcing  process 
(Berloff,  2005).  Therefore,  the  results  obtained  in  the  present 
study  may  be  useful  for  interpretation  of  a  wide  spectrum  of 
problems  related  to  model  predictability  in  the  atmosphere 
and  ocean. 

The  rest  of  the  paper  is  organized  as  follows.  Section  2 
explains  a  predictability  metrics  used  to  quantify  the  model 
predictability  for  both  small-  and  large-amplitude  perturba¬ 
tions.  Features  of  the  reference  solution  (the  control  run)  are 
discussed  in  Sect.  3.  The  surface  wind  is  decomposed  into 
two  parts:  steady  (“climatic”)  part  and  stochastic  one  caused 
by  unknown  synoptic  variability.  Statistics  of  the  wind  un¬ 
certainty  is  given  in  Sect.  4.  The  model  phase  space  is  intro¬ 


duced  in  Sect.  5.  Section  6  investigates  the  sensitivity  of  the 
root  mean  square  prediction  error  to  variations  of  stochas¬ 
tic  wind  and  the  tolerance  level  (the  accepted  prediction  ac¬ 
curacy).  Section  7  analyzes  finite-amplitude  induced  phase 
transitions  of  predictability.  Section  8  develops  a  technique 
for  identification  of  the  probability  density  function  of  the  ir¬ 
reversible  predictability  time.  Herein,  we  demonstrate  that 
the  statistics  of  this  time  is  rather  Weibullian  than  Gaussian. 
The  predictability  horizon  is  estimated  in  Sect.  9.  Section  10 
provides  the  conclusions.  Appendix  A  contains  analytical 
representations  for  the  wind  error  source  term. 


2  Predictability  measures 


We  examine  the  sensitivity  of  a  reference  solution  relative 
to  stochastic  variations  of  wind  forcing.  Such  sensitiv¬ 
ity  may  be  measured  by  the  traditional  non-dimension  root 
mean  square  difference  between  perturbed  (“pert”)  and  non- 
perturbed  (“ref”)  solution  presented  by  a  variable  T,  which 
may  stand  for  energy,  temperature,  salinity,  the  stream  func¬ 
tion  etc.  (Holland  and  Malanotte-Rizzoli,  1989;  Brasseur 
et  al.,  1996;  Robinson  et  al.,  1996;  Wirth  and  Ghil,  2000, 
among  others), 

<  7(f)  >=<  /2ert(f)  >  /4f(f)  (1) 

and  the  irreversible  predictability  time  (IPT)  (Ivanov  et  al., 
1994)  defined  as 


r  (e)  =  inf 

t>  0 


>  e2 


)• 


(2) 


where  /pert=  || ^pert — ^ref || ,  /ref=  ll^refll,  hereafter  <...>  is 
the  ensemble  averaging,  s  —  e//ref  is  the  non-dimension 
tolerance  (the  accepted  prediction  accuracy),  ||  ||  is  the  Eu¬ 
clidean  norm.  According  to  Eq.  (2)  the  irreversible  pre¬ 
dictability  time  is  a  time  at  which  the  prediction  error  71-/2 
reaches  a  predetermined  level  e  for  the  first  time,  i.e.  any  re¬ 
turns  of  model  predictability  are  impossible  after  r(e). 

The  IPT  is  clearly  interpretable  in  a  model  phase  space, 
where  perturbed  and  non-perturbed  (reference)  solutions  are 
represented  by  X'  and  X  trajectories,  respectively,  and  the 
equation  I{t)—sr  describes  a  spheroidal  surface  S(e)  mov¬ 
ing  along  the  trajectory  X  (Fig.  la).  A  distance  \Xa— X'a\ 
between  the  trajectories  at  time  moment  tQ  (the  initial  er¬ 
ror)  usually  grows  with  time  due  to  model  inaccuracy,  and 
becomes  larger  than  s  (crossing  S(s))  after  a  time  r(e) 
(Fig.  la).  This  time  is  defined  as  the  IPT. 

For  a  steady  reference  solution  (represented  by  a  fixed 
point  in  the  model  phase  space),  the  IPT  becomes  the  clas¬ 
sical  first  passage  time  (FPT).  The  classical  FPT  is  a  time  at 
which  a  trajectory  reaches  a  boundary  for  the  first  time  (Gar¬ 
diner,  1985).  Therefore,  the  IPT  can  also  be  defined  as  the 
FPT  for  varying  boundaries  (compare  Fig.  la  and  Fig.  lb). 
The  FPT  plays  an  essential  role  in  many  applied  fields.  We 
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Fig.  1.  Definitions  of  (a)  IPT  and  (b)  FTP.  Phase  trajectories  X 
(the  reference  solution)  with  the  initial  position  of  Ag,  and  X'  (a 
perturbed  solution)  with  the  initial  position  of  X^are  denoted  by 
solid  and  dashed  curves,  respectively.  A  time  when  X'  crosses  5(e) 
in  the  point  A  is  the  IPT  (a)  or  FPT  (b). 

can  suppose  that  the  IPT  may  become  a  useful  tool  for  the 
analysis  of  ocean  model  predictability  too. 

Statistics  of  the  IPT  (Ivanov  et  al.,  1994;  Chu  and  Ivanov, 
2005)  can  be  represented  by  the  probability  density  func¬ 
tion  (r-PDF)  or  the  cumulative  distribution  function  (r-CDF) 
[P(e,  f— to)]-  r-CDF  is  the  probability  that  r>t—to  for  a 
given  tolerance  level. 

In  practical  applications,  statistics  of  IPT  (r-CDF  or  r- 
PDF)  may  be  identified  from  r-moments, 

00 

r /  (e)  =  /  J  (t  —  ?o)/_1  P  (e,  t  —  to)dt,  1=1,...,  L.  (3) 

to 

Knowing  these  moments  we  may  compute  r-mean  ,  r- 
variance  (r-var)  ,  r-skewness  (r-sk)  and  r-kurtosis  (r-ku). 
For  example, 

r-mean  =  n,  r-var  =  t2  —  (ri)2  etc. 

If  the  cumulative  distribution  function  P  has  a  heavy  tail 
for  large  values  of  t— to,  high-order  r-moments  (1—2,  ...,  L) 
sometimes  do  not  exist  because  the  integral  in  Eq.  (3)  does 
not  converge.  In  this  case  we  suggested  (Ivanov  and  Chu, 
2007)  to  identify  r-CDF  or  r-PDF  from  the  probability- 
weighted  moments  ( a / )  defined  originally  by  Greenwood  et 
al.  (1979), 

1 

oil  —  J  X(P)  (\  —  P)]  dP ,  l  =  l,...,L  ,  (4) 

o 

where  X(P)  is  the  quantile  function  ( i.e.,  the  inverse  of  cu¬ 
mulative  distribution  function). 


Fig.  2.  Basin  geometry.  The  x\  and  %2  axes  point  toward  east  and 
north,  respectively. 


In  practice  the  same  moments  aj  are  estimated  from  an  or¬ 
dered  random  sample  (r)^=1  of  size  N  (Hosking  and  Wallis, 
1997) by 

n=  1 

where  C"_/  are  the  binomial  coefficients. 

The  probability-weighted  moments  always  exist  and  are 
robust  relative  to  sampling  error.  Therefore,  the  robust  es¬ 
timate  of  r-CDF  or  r-PDF  for  small  forecast  ensembles  is 
possible.  This  is  one  of  advantages  of  the  IPT.  An  appropri¬ 
ate  method  for  estimating  distribution  functions  from  knowl¬ 
edge  of  the  probability  weighted  moments  will  be  discussed 
in  Sect.  8. 

The  moments  of  IPT  as  functions  of  X0—X'0  satisfy  the 
Pontryagin-Kolmogorov-Stratonovich  equations  (Pontryagin 
et  al.,  1969),  which  are  linear  elliptic  differential  equations. 
Their  asymptotic  solutions  can  be  obtained  in  many  cases. 
For  example,  Chu  et  al.  (2002)  calculated  analytically  first 
two  moments  of  IPT  for  a  low-order  nonlinear  atmospheric 
model  (Lorenz,  1984).  Therefore,  the  analytical  estimate  of 
model  predictability  in  the  IPT  context  is  another  advantage 
of  our  approach. 


3  The  reference  solution 

We  consider  a  rectangular  semi-closed  basin  with  the  hori¬ 
zontal  dimensions:  L  i = 1  050  k m  and  7-2=1000 km,  and  with 
constant  depth  H= 2  km,  which  is  situated  on  a  mid-latitude 
/1-plane.  The  basin  has  rigid  (T)  and  open  (T')  boundaries. 
The  geometry  of  the  basin  and  its  sizes  are  shown  in  Fig.  2. 
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Our  numerical  model  is  the  nonlinear  shallow-water  equa¬ 
tions  with  nonlinear  bottom  friction,  wind  and  boundary 
forcing 

d  Du  |  \  it  , 

- - +L{Du2 ,  Dui)-fDu2=-gDW^+Wl-aEl/2ul,  (6) 

dt 

~^+L(Dui,  Du2)+fDul=-gDV2t;+W2-aE1l2U2 ,  (7) 


and  the  mass  conservation  equation 


dt 

— — h  (V[Du\  +  V2Z)m2)  —  0, 
at 


(8) 


where  L(...,  ...)  is  the  nonlinear  advective  operator; 
[Vi,V2]=[gfr,  ] ;  m  1  and  u2  are  the  zonal  and  meridional 
velocities,  respectively;  D=H+t;,  t;  is  the  sea  surface  ele¬ 
vation;  the  drag  coefficient  cn=2.5xl0-3;  the  gravity  g:  and 
E—u\+u2. 

The  Coriolis  parameter  varies  linearly  with  a  beta 
plane  approximation  /=/o+0x2,  where  fy=2Q  sin (cp0)  and 
P=(2£l/a)  cos((p0).  Here,  £2  and  a  are  the  rate  of  ro¬ 
tation  and  the  radius  of  the  Earth,  respectively;  ^o=35°. 
For  the  chosen  model  parameters:  /O=7.3xl0-5  s-1, 
0=2. Ox  10-11  m-1  s-1. 

A  flow  in  the  semi-closed  basin  bounded  by  T  U  T'  is 
forced  by  both  the  zonal  wind  forcing  W 1  (W2—0)  varying 
with  latitude  as 


W\  = 


ws  ttx  2 

- cos( - ), 

pw  L2 


(9) 


where  /?u,=1025  kgm-3,  ws  is  the  wind  stress, 
^L=1.0x  10-3  m2  s-2  .and  a  prescribed  net  flux  (char¬ 
acterized  by  the  normal  velocity  ui,(x2.  t)  and  surface 
elevation  gb{x2,t)  along  the  boundary  E').  Zero  normal 
velocity  and  zero  Neumann  conditions  for  the  surface 
elevation  are  imposed  on  the  rigid  boundary  T. 

The  chosen  model  configuration  is  suitable  for  the  analysis 
of  ocean  model  predictability  affected  by  different  kinds  of 
stochastic  uncertainties:  errors  inserted  in  initial  conditions 
(Ivanov  and  Chu,  2007),  wind  (the  present  study)  and  open 
boundary  conditions1.  Cross-correlations  between  these  er¬ 
rors  can  also  be  studied. 

Model  (6-8)  is  similar  to  that  used  by  Veronis  (1966)  for 
the  analysis  of  nonlinear  wind-driven  circulation  in  a  closed 
basin.  But  in  contrast  to  Veronis  (1966)  we  parameterize  bot¬ 
tom  friction  by  the  quadratic  drag  law  (Pedlosky,  1987). 

The  barotropic  mode  of  Princeton  Oceanographic  Model 
(Blumberg  and  Mellor,  1987)  was  applied  to  Eqs.  (6-8)  with 
the  following  model  parameters:  spatial  resolution  -  50  km; 
time  step  -  2  min. 

The  prescribed  non-stationary  net  flux  across  the  open 
boundary  is  computed  as  it  was  explained  in  Chu  et 

1  Ivanov,  L.  M.  and  Chu,  P.  C.:  Effects  of  stochastic  open  bound¬ 
ary  uncertainty  on  predictability  of  regional  ocean  models,  Mon. 
Weather  Rev.,  in  preparation  ^^Hstatus?,  2007.) 


al.  (1997).  The  structure  of  open  boundary  conditions  on 
day-0  and  day-60  is  demonstrated  in  Figs.  3a  and  b,  respec¬ 
tively.  The  initial  condition  represents  a  non-closed  anti- 
cyclonic  gyre  shown  in  Fig.  3a.  The  corresponding  initial 
surface  elevation  is  not  shown  because  its  structure  is  obvi¬ 
ous. 

After  30  days  of  integration  the  model  reaches  a  spin  up 
when  the  spatially  averaged  kinetic  energy  oscillates  with  a 
period  of  120  days.  Amplitude  of  this  oscillation  reduces 
with  time  exponentially  with  rate  of  1000  day-1.  The  spa¬ 
tially  averaged  kinetic  energy  for  the  first  60  days  is  shown 
in  Fig.  4a  only  because  this  time  period  is  used  for  sensi¬ 
tivity  studies.  The  circulation  pattern  formed  after  day-30 
presents  a  multi-gyre  structure  with  maximum  velocities  up 
to  0.9..1.0m/s  (Fig.  3b)  and  high  surface  elevation  near  lm 
(not  shown). 


4  Wind  forcing  uncertainty 


Governing  Eqs.  (6-8)  are  perturbed  by  adding  the  stochastic 
wind  forcing  w—  (uq,  w2)  to  W—  (W\,  W2).  The  stochastic 
wind  forcing  is  traditionally  parameterized  in  the  following 
form 

w(x  1,  xt,  t)  =  ^-Cd  | C/(:ri,  xt,  t) |  U (xi,  xt,  t),  (10) 

Pw 

where  pajr  is  the  air  density  (1.3  kgm-3),  Cj(2  x  10-3)  is  the 
drag  coefficient,  and  U is  the  stochastic  wind. 

Following  Sura  et  al.  (2001)  U  is  represented  by 

V—  [U\{x\,  X2,  t),  U2(x  1,  x2,  ?),]  =/r(f)erG1/2(xi,  x2),  (11) 


where  [p,\ (t),  p2(t)]  are  white  Gaussian  vector  pro¬ 

cesses  with  zero  mean  and  unit  variance;  a2  is  the  wind  vari¬ 
ance;  the  spatial  structure  function  G  characterizes  a  degree 
of  spatial  inhomogeneity  of  wind  perturbations  above  an  area 
of  interest. 

Two  different  structure  functions  G  are  used.  The  first  one 
is  given  by 


Gi{xi,xt_)  =  cos(— -).  (12) 

E2 

In  this  case  only  the  amplitude  of  wind  stress  (Eq.  9)  is  dis¬ 
torted  by  the  non-Gaussian  white  noise. 

The  second  one  is  chosen  as 


G2(x  1,  X2)  =  Q! scale 


^ifterfl  -r  I  erf  I 


20i 


202 


-1/2 


exp  -- 


(xi-Li/2)2  (x2-LtJ2Y 


202 


202 


■  (13) 


Here,  erf  is  the  error  function;  ascaie  is  a  scaling  parame¬ 
ter;  (0i,  02)  are  the  decorrelation  scales;  G2  shows  the  im¬ 
pact  of  the  localized  atmospheric  eddy  activity  near  the  point 
(L i/2,  L2/ 2)  on  the  surface  wind  perturbations  (Sura  et  al., 
2001). 
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Fig.  3.  Spatial  structure  of  the  reference  solution  at  the  initial  state  (a)  and  after  integration  for  60  days  (b).  Open  boundary  conditions  (u/,) 
corresponding  to  the  reference  solution  are  shown  to  the  right  of  the  circulation  patterns. 


In  most  numerical  experiments  the  scaling  constant  ascaie 
is  chosen  to  adjust  the  weight  function  in  Eq.  (13)  to  1  for 
/3c.=y81=/32=600km.  However,  a  number  of  computations 
use  pc  between  100  km  and  600  km. 

The  noise  in  the  surface  wind  with  cr2=28.0m2  s-2  cor¬ 
responds  to  typical  observed  atmospheric  conditions  in  the 
North  Atlantic  region  (Wright,  1988).  Therefore,  the 
stochastic  forcing  (Eqs.  10-13)  is  a  conceptual  tool  to  study 
the  effect  of  noise  on  simple  and  more  complex  wind-driven 
regional  ocean  models. 

To  understand  statistics  of  wind  forcing,  Eq.  (10)  is  re¬ 
written  into 

w(xi,x2,  t)  =  —  Q<t2G(xi,x2)  \fi(t)\n(t) 

Pw 

=  —  Cdo2G{xi,x2)w,  (14) 

Pw 

where  w—  pAt).  Then,  the  probability  density  function 
f(w i,  wj)  is  calculated  using  the  elementary  zero-memory 
transformations,  which  are  discussed  in  most  textbooks  of 
probability  theory  (see,  for  example,  Stratonovich,  1963). 
Accordingly  to  the  general  theory 

f(u> t,  u>2 )  =  /[gf'ODi,  w2),  g2 u>2)]  ■  |7| ,  (15) 


where  J  is  the  Jacobian  of  the  transformation  from  the  ran¬ 
dom  variables  // 1  and  /ii  to  the  random  variables  ih\  and  w2\ 
gj-1  and  gj~  1  are  the  inverse  functions. 

Simple  calculations  result  into 


f(wi,  w2)= 


4: r  (wj+it^) 


2x1/2 


exp 


-  (w\+w^  /  /2 


■(16) 


The  means  (wi),(u>2)  and  variances  cr2,cr2  computed  from 
Eq.  (16)  have  the  following  values 


(w\ )  =  (w2)  —  0.0  and  =a\  =  3.0.  (17) 


Both  and  w2(t)  are  delta-correlated  processes  (Kly- 

atskin,  2005). 


Fig.  4.  Characteristics  of  the  reference  solution:  (a)  the  kinetic 
energy  averaged  over  the  semi-closed  basin,  (b)  the  relative  variance 
Sm  computed  for  the  initial  state  of  the  reference  solution  (solid 
curve)  and  after  integration  for  60  days  (dashed  curve). 


In  the  polar  coordinate  system  {z,  0}[wi=s  cos(0), 
w2=z  sin($)],  probability  density  function  (Eq.  16)  trans¬ 
forms  to  the  following  form: 

f(z,  61)  =  7-  exp  (-z/2) .  (18) 

47T 

Two-dimensional  distribution  Eq.  (18)  is  easily  transformed 
to  a  one -dimension  form  /;nt  simply  by  integrating  f(z,  0) 
with  respect  to  0 , 

/int(z)  =  ~  exp(  z/2),  (19) 

which  is  the  exponential  distribution. 

The  above  calculations  clearly  indicate  that  although  U  is 
a  white  Gaussian  process,  the  wind  stress  w  is  not.  Therefore 
PE  has  non-Gaussian  statistics  even  if  the  wind  uncertainty  is 
small.  For  large  values  of  stochastic  wind  stress,  distribution 
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Fig.  5.  Four  orthonormal  modes  i jrm  with  numbers  m= 1,3,9,  and  30.  The  contour  interval  is  non-dimensional,  with  positive  vorticity  in  dark 
and  negative  vorticity  in  light,  and  non-dimensional  velocity  vectors  are  overlaid  in  each  panel. 


function  (19)  decays  slower  than  the  Gaussian  one.  This  indi¬ 
cates  that  rare  high-energetic  wind  events  can  strongly  con¬ 
tribute  to  w  statistics.  Therefore,  even  for  small-amplitude 
errors,  r-PDF  has  asymmetric  shape  with  a  tail  stretching 
into  short  prediction  time  scales. 

A  weak  second-order  algorithm  (Cao  and  Pope,  2003) 
is  applied  to  numerical  integration  of  stochastically  forced 
Eqs.  (6-8).  Although  the  time  step  for  the  model  integra¬ 
tion  is  5  min,  the  stochastic  wind  is  updated  every  hour.  A 
correlation  time  (tc)  for  such  a  noise  is  shorter  than  one 
hour.  Since  characteristic  time  scales  for  the  reference  so¬ 
lution  and  PE  are  rref~  1 0-20  days  (determined  from  Fig.  4a) 
and  terror^ 3—5  days  (determined  from  Fig.  7b),  respectively, 
rmtl  k / terror «1  •  Such  a  stochastic  wind  represents 
a  white  noise-like  process  (Stratonovich,  1963). 

PE  statistics  is  insensitive  to  more  often  update  of  the 
stochastic  wind.  We  made  computations  with  update  varying 
from  less  than  one  hour  to  5  min.  These  computations  re¬ 
quired  very  large  computer  resources.  Therefore,  the  choice 
of  one-hour  update  is  a  trade-off  between  accuracy  in  repre¬ 
sentation  of  the  wind  forcing  and  the  computational  cost. 

Ensembles  of  perturbed  model  trajectories  were  used  to 
compute  r-PDF.  Little  difference  in  r-statistics  is  obtained 


between  ensembles  of  1  x  103,  5xl03,  lxlO4,  2xl04,  and 
5xl04  samples.  The  optimal  size  of  an  ensemble  sampling, 
i.e.  a  number  of  ensemble  realizations  providing  a  trade¬ 
off  between  the  ensemble  ability  to  reproduce  main  features 
of  PE  statistics,  and  the  computational  cost,  is  estimated 
as  103  for  any  values  of  er2.  Flereafter  the  non-dimension 
variance  of  wind  perturbations  introduced  as  a2—o2/c rj, 
a(j—  1 .0  m2  s~2,  is  used.  The  optimal  size  is  found  using  the 
non-symmetrical  Kullback-  Leibler  distance  (White,  1994). 
An  interested  reader  is  referred  to  Chu  and  Ivanov  (2005), 
Ivanov  and  Chu  (2007)  for  more  details  of  such  an  approach. 

5  Model  phase  space 

For  the  chosen  model  parameters,  a  quasi-geostrophic  ap¬ 
proximation  (Pedlosky,  1987)  is  applicable  to  interpret  the 
reference  and  perturbed  flows  in  a  model  phase  space.  The 
basis  of  M-dimension  phase  space  is  formed  from  orthonor¬ 
mal  functions  (modes)  \j/m,  which  are  the  eigenvectors  of  the 
plane  Laplace  operator  V2  (Eremeev  et  al.,  1991,  1992), 

^  /km  =  t/bn  IrUT'  =  0,  Ilf  =  1,  ...,  M.  (20) 
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Then,  the  geostrophic  stream  function  is  decomposed  as 
M 

f(X\,X2,  t)=  £  Am(t)\lrm(x l,X2)+ifhmn(xi,X2,  t)+C(t), 

m= 1 

M  =  100  ,  (21) 

where  i/Tarm  is  the  harmonic  function;  the  constant  C(t)  is 
determined  from  the  mass  conservation  constraint  imposed 
upon  the  stream  function  (McWilliams,  1977). 

The  spatial  modes  \j/m  are  determined  only  by  the  geom¬ 
etry  of  a  basin,  and  can  be  easily  computed  for  any  non- 
rectangular  domain.  Figure  5  shows  the  spatial  structure  of 
a  few  basis  functions  i jrm  involved  in  the  present  analysis. 
The  spatial  modes,  in  general,  have  no  physical  significance 
by  themselves,  but  only  when  they  imply  a  flow.  However, 
they  are  useful  to  identify  the  energy-dominated  scales  for 
the  reference  and  perturbed  flows. 

The  harmonic  function  V^harm  is  obtained  from 

V^l/'liarm  =  0  ,  l/^harm  |r  —  0  , 

x 2 

lAharm  lr'  =  -  J  ub(t ,  y)dy ,  (22) 

0 

It  is  a  highly  predictable  component  of  the  flow  because  of 
the  exact  value  of  ii],  in  Eq.  (22). 

The  PE  presents  the  sum  of  the  mean  or  systematic  error 
Vhef—  (V^pert)  and  the  transient  or  random  error  8 x[r: 

[ipertj  ~  (  II  Vhef  —  (V^pert)  —  W  ||  j 

=  ||^ref-(^pert)||2+(ll^||2).  (23) 

The  inertial  (nonlinear)  terms  of  the  governing  equations 
hardly  contribute  to  ( V^pert)  at  the  initial  stage  of  PE  growth, 
and  their  contribution  is  negligible  at  later  stages.  The  ran¬ 
dom  error  grows  faster  than  the  systematic  error.  Therefore, 
we  suggest  quantifying  the  PE  behavior  through  the  growth 
of  the  random  error  only. 

The  reference  and  perturbed  solutions  are  repre¬ 
sented  in  the  model  phase  space  as  the  reference 
A=  Am(0]  and  error  a—  [ai(f), ...,  om(0] 

trajectories,  respectively.  Using  these  notations  the  variance 
of  the  random  prediction  error  becomes 

<ll«*ll>  =  £(4).  (24) 

m=  1 

and  the  wind  error  source  term  (see  Appendix  A)  is  written 
by 

Rm  =  Ym  (bl  +  4)  ,  (25) 

where  ym=^CdH~lk~la2,  bm=  ff  xlrmdx\dx2  , 
cm—  ff  i/mdx\dx2,  the  double  integration  is  made  over 
the  semi-closed  basin  area. 
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Fig.  6.  Phase  portrait  in  a  phase  sub-space  generated  by  the  basis 
functions  xjfq  and  i //s. 


Applying  the  classical  linear  stability  analysis  (Gucken- 
heimer  and  Holmes,  1983)  to  the  model  (6-8)  linearized 
near  the  spin  up  solution,  we  find  that  the  spin  up  represents 
an  unstable  focus  (spiral  point)  in  the  model  phase  space. 
The  growth  and  decay  of  infinitesimal  perturbations  near  this 
point  are  characterized  by  the  spectra  of  positive  and  negative 
local  characteristic  exponents  only. 

Therefore,  an  error  trajectory  should  tend  to  this  focus  (de¬ 
noted  by  B  in  Fig.  6)  along  stable  manifolds  corresponding  to 
the  large-scale  negative  exponents  and  simultaneously  drifts 
from  B  along  the  small-scale  unstable  manifolds  correspond¬ 
ing  to  the  positive  exponents.  Such  a  model  trajectory  is 
asymptotically  unstable  in  Lyapunov  sense  as  t—>o o  (Guck- 
enheimer  and  Holmes,  1983). 

Figure  6  shows  projections  of  error  trajectories  onto  the 
phase  subspace.  The  trajectories  tend  to  reach  the  focus 
along  the  stable  manifolds  projected  onto  the  phase  plane 
[cip,  ciq].  However,  they  move  away  from  the  focus  along  un¬ 
stable  manifolds  projected  on  the  phase  planes  [ap,  av]  and 
[aq,  as ]. 

We  use  M=100  and  confirm  that  such  a  choice  does  not 
smooth  the  reference  trajectory  for  70  days  of  model  integra¬ 
tion.  The  relative  variance  Sm=  {Ap-Ap)mp_\ '  Iff  converges 
to  1  very  quickly  as  m  increases  (Fig.  4b).  The  first  fifty-sixty 
modes  contain  more  than  99%  of  variance  for  the  reference 
solution. 

One  hundred  mode  representation  is  also  quite  suf¬ 
ficient  to  approximate  the  error  trajectory  for  60- 
70  days  of  model  integration.  The  relative  variance 

S'm=  •  lylpeitj  converges  to  1  as  m  in¬ 

creases,  slower  than  Sm  (compare  Figs.  4a  and  7a),  but  the 
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Fig.  7.  Features  of  a  perturbed  solution:  (a)  the  relative  variance 
S'm  computed  at  day-10  (dashed  curve)  and  day-50  (solid  curve), 
cr2=1.0;  (b)  </pert>  normalized  by  its  saturation  value  </2 > 
for  cr2=0.1  (solid  curve),  1.0  (dotted  and  dashed  curve),  and  2.0 
(dashed  curve). 


Fig.  8.  The  growth  rate  Q  (solid  curve)  for  different  a~  and  G:  (a) 
a2  =  0.1,  G=G2;  (b)  o-2=1.0,  G=Gi;  (c)  a2=1.0,  G=G2,  and 
(d)  cr2=2.0,  G=G2.  Dashed  line,  white  dots  and  asterisks  show 
exponential,  linear,  and  power  (with  scaling  exponent  of  8.8  x  10“ 1 ) 
laws,  respectively. 


speed  of  the  convergence  is  quite  high. 


6  Error  evolution 

Typical  growth  of  an  ensemble  averaged  PE  with  time  for 
G—G i  and  cr2  varying  between  0.1  and  2.0,  is  given  in 
Fig.  7b.  Perturbations  excited  by  uncertainty  of  stochastic 
wind  grow  at  all  scales  and  during  the  whole  50-60  day 
period.  That  is  in  contrast  to  the  case  when  there  is  un¬ 
certainty  in  the  initial  condition  only.  In  the  last  case  the 
high  predictability  of  the  dynamical  regime  within  the  initial 
15-20  day  period  was  clearly  demonstrated  by  Ivanov  and 
Chu  (2007):  the  PE  at  first  decays  with  time  for  all  scales 
due  to  dissipation  caused  by  nonlinear  bottom  friction,  and 
only  after  day -20  grows  faster  than  [quasi] -exponentially. 
Therefore,  the  presence  of  the  spatio-temporal  noise  (Eq.  10) 
in  wind  forcing  (Eq.  9)  causes  the  monotonic  error  growth 
shown  in  Fig.  7b. 

At  least  four  predictability  regimes  are  identified  from 
Fig.  7b.  In  all  the  cases  the  PE  grows  in  a  monotonic 
manner  but  with  different  speeds.  More  accurately,  these 
regimes  can  be  identified  using  the  growth  rate  defined  as 

2=37ln</pert>- 

A  set  of  growth  rates  computed  for  different  G  and  cr2 
is  presented  in  Figs.  8a,  b,  c,  d.  These  results  clearly  indi¬ 
cate  that  error  dynamics  strongly  depends  on  the  intensity 
and  spatial  inhomogeneity  of  wind  uncertainty. 


6.1  Linear  growth  of  perturbations 

At  the  initial  stage  (transient  phase)  where  the  stochastic 
forcing  term  dominates  the  governing  equations,  Q~\/t 
(Figs.  8a,  b,  c,  d).  It  corresponds  to  the  linear  growth  of 
the  mean  square  error: 

(/pert)  *  DeSt,  (26) 

Duration  of  this  regime  is  typically  up  to  4-5  days  if  d2~0.1- 
1.0.  The  effective  coefficient  De g-  is  determined  by  sum¬ 
mation  of  contributions  from  the  error  source  term  at  all 
M 

wavenumbers  /7en—  Rm- 

m=  1 

Linear  law  Eq.  (27)  was  earlier  documented  in  a  num¬ 
ber  of  studies  (see,  for  example  Vannitsem  and  Toth,  2002). 
We  have  analytically  determined  the  wind  error  sources  for 
model  (6-8)  (Appendix  A).  Our  calculations  show  strong  de¬ 
pendence  of  the  effective  coefficient  on  the  variance  of  wind, 
as  ~cr4,  as  well  as  on  degree  of  spatial  inhomogeneity  of  the 
wind  forcing. 

6.2  Power  growth  of  perturbations 

For  moderate  but  inhomogeneous  winds  the  power  growth  of 
perturbations  are  observed  in  our  numerical  experiments  af¬ 
ter  the  transient  phase  (for  example,  see  Figs.  8a,  b,  c).  For 
small  values  of  (t2<$(L0  the  power  growth  is  replaced  by  the 
exponential  growth  (shown  by  the  dashed  line  in  Fig.  8a). 
If  cr 2  exceeds  1,  there  is  no  exponential  growth  and  the  PE 
grows  with  the  power  law  with  power  exponent  of  about 
8.8  xl0_1.  This  regime  exists  between  day-5  and  day-15  in 
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Fig.  9.  Spectra  of  wind  error  term  with  (a)  G j,  (b)  Gj  and 
A=600  km.  and  (c)  Gj  and  X=  100  km.  A  critical  wavenumber  is  de¬ 
termined  from  the  condition  of  Rm/  max(/?,H)=1.0x  10-2  (shown 
by  dashed  line).  Black  arrows  indicate  the  critical  numbers. 


Fig.  10.  Coherent  behavior  of  perturbations  on  different  scales  for 
ff2=1.0  and  G=G2-  (a)  Large,  and  (b)  small  scale  perturbations. 
The  largest  scale  perturbation  with  m= 1  is  labeled  by  dashed  curve. 


Fig.  8a,  between  day-7  and  day-23  in  Fig.  8b,  between  day- 
4  and  day- 14  in  Fig.  8c,  but  there  is  no  power-law  regime 
in  Fig.  8d  when  the  stochastic  wind  uncertainty  is  too  large 
("o'2  >2.0).  For  such  a  variance  the  linear  growth  of  PE  dom¬ 
inates. 

For  G—G  i  the  spectrum  of  Rm  is  linear  with  dominat¬ 
ing  peaks  at  wavenumbers  m=l,3,  and  5  (Fig.  9a).  These 
wavenumbers  indicate  modes  with  maximum  response  to  the 
stochastic  wind  forcing.  The  weak  wind  forcing  (er2<^1.0) 
essentially  affects  the  large  scales  of  the  flow  and  excites 
only  several  low-order  modes.  In  this  case  the  PE  at  first 
grows  linearly,  then  its  quasi-exponential  growth  is  observed. 
Smaller  scales  affected  by  the  stochastic  wind  are  subject 
to  strong  viscous  damping  due  to  increasing  drag  coefficient 
a  with  growth  of  the  kinetic  energy  of  large-scale  perturba¬ 
tions.  Therefore,  the  smaller  scales  grow  slower  than  the  un¬ 
stable  large  scales.  The  growing  perturbations  rapidly  adopt 
the  horizontal  scales  comparable  to  those  of  the  reference 
state. 

Alternatively,  stronger  stochastic  wind  (a 2  >1.0)  excites 
more  modes  at  smaller  scales  than  the  weak  wind,  and  the 
coherent  behavior  of  modes  is  clearly  observed  in  this  case 
(for  example  see  Fig.  8b). 

For  G—G 2  and  /Jc=600km,  what  corresponds  to  inho¬ 
mogeneous  winds,  the  spectrum  of  Rm  is  continuous  and 
band  limited  at  the  critical  wavenumber  ;;;=11  (Fig.  9b). 
Therefore,  the  stochastic  wind  excites  several  modes  around 
wavenumbers  of  2  and  10.  The  PE  growth  ratio  depends  on 
the  spectrum  of  local  characteristic  exponents.  In  this  case 
the  coherent  behavior  of  modes  exists  even  for  the  weak  wind 
forcing  (ct2<$(1.0). 

The  decorrelation  scale  fic  determines  the  width  of  spec¬ 
trum  band  for  Rm.  For  example,  the  critical  wavenumber 


equals  to  32  if  ftc=  1 00  km  (Fig.  9c).  Reduced  values  of  fic 
lead  to  a  wider  spectrum  of  local  characteristic  exponents 
and  stronger  contribution  of  the  cumulative  effects  to  the  PE 
growth. 

Our  computations  have  also  shown  that  one  of  the  main 
specificities  of  the  power  growth  regime  is  that  all  dominat¬ 
ing  scales  (modes)  may  exhibit  a  similar  growth  rate.  For 
example.  Fig.  10  demonstrates  two  groups  of  modes  with  dif¬ 
ferent  coherent  behavior  from  day -7  to  day-25.  The  coherent 
behavior  of  modes  is  a  collective  response  to  the  external 
stochastic  forcing  when  amplitudes  of  many  modes  exceed 
some  threshold  at  the  same  time  due  to  spatial  inhomogene¬ 
ity  of  wind  forcing  (defined  by  the  choice  of  G).  Clearly  that 
the  coherent  behavior  is  absent  if  this  threshold  was  exceeded 
only  by  a  few  modes. 

6.3  “Super-exponential”  growth  of  perturbations 

After  day-25  the  PE  grows  faster  than  exponentially  (“super- 
exponentially”)  until  non-linear  interactions  between  differ¬ 
ent  scales  destroy  this  growth  (saturation  regime)  (Fig.  7b). 
Strong  coherence  in  behavior  of  different  modes  accompa¬ 
nies  the  “super-exponential”  growth  of  perturbations. 

Figure  11a  shows  coherent  behavior  of  30  dominated 
modes.  Explicit  “synchronization”  in  behavior  of  these 
modes  is  observed.  After  35-37  days  of  integration,  the  first 
mode  (shown  by  dashed  curve  in  Fig.  10a)  is  a  “driver”  deter¬ 
mining  the  behavior  of  all  other  dominant  large-scale  modes, 
which  are  called  “responses”  (hereafter  we  use  the  terminol¬ 
ogy  from  Boccaletti  et  ah,  2002).  During  a  10-12  day  time 
period  (up  to  day-47)  the  driver  and  responses  perform  coher¬ 
ently.  For  small-scale  perturbations  at  least  two  drivers  can 
be  identified  in  Fig.  lib.  Mode  m= 9  (the  first  driver)  grows 
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Fig.  11.  Coherent  behavior  of  perturbations  on  different  scales  for 
cr2=0.1  and  G=G2-  (a)  Mode  m= 1  is  “the  driver”  for  large-scale 
responses.  Three  responses  (modes  m=  4,  12,  29)  are  indicated  by 
bold  dashed  lines;  (b)  two  drivers  (modes  m= 9,  11)  observed  for 
small-scale  perturbations.  Mode  m= 9  grows  with  the  exponential 
law. 


along  the  exponential  law.  The  non-exponential  growth  of 
the  second  driver  (mode  »t=ll)  is  a  consequence  of  the 
fact  that  different  scales  grow  with  different  local  char¬ 
acteristic  exponents  varying  between  2.3xl0-1  day-1  and 
9.4xl0-1  day-1. 

For  very  strong  external  noises  the  modes  may  not  be  syn¬ 
chronized  because  the  intensive  noise  destroys  correlations 
among  them.  This  effect  is  clearly  observed  in  our  numerical 
experiments  independently  on  G  if  rf 2  >2.0. 

6.4  Physical  mechanisms  of  perturbation  growth 

The  numerical  results  discussed  above  explicitly  show  “a 
synchronization”  effect  due  to  external  noise  and  replace¬ 
ments  of  the  traditional  exponential  growth  of  perturbations 
by  the  power  or  sub-exponential  growth.  There  is  a  num¬ 
ber  of  mechanisms  for  which  noise  can  lead  to  more  order 
in  the  dynamics.  To  be  mentioned  here  are  the  effects  of 
noise-induced  order  in  chaotic  dynamics  (Matsumoto  and 
Tsuda,  1983),  synchronization  of  self-sustained  oscillators 
(Pikovsky  et  al.,  2000),  cumulative  effects  of  many  different 
scales  (Aurell  et  al.,  1996),  coherence  resonance  (Pikovsky 
and  Kurths,  1997),  stochastic  resonance  (Nicolis  and  Nico- 
lis,  1981;  Benzi  et  al.,  1981),  and  interference  between  initial 
error  and  stochastic  forcing  (Seki  et  al.,  1993).  These  effects 
in  some  respects  are  close  and  cannot  be  easily  distinguished 
from  one  another  when  signals  reflect  different  variables  of 
the  same  system  (Rosenblum  et  al.,  2004). 

In  our  case  circulation  dynamics  is  not  driven  by  a  peri¬ 
odical  force.  That  allows  excluding  the  stochastic  resonance 
as  a  possible  physical  mechanism  driving  mode  dynamics. 


The  stochastic  resonance  appears  if  both  periodic  and  noisy 
forces  drive  a  nonlinear  system,  with  the  periodic  response 
having  a  maximum  at  some  noise  amplitude. 

Perez -Munuzuri  et  al.  (2005)  have  demonstrated  a  coher¬ 
ent  resonant  behavior  for  an  atmospheric  global  circulation 
model  induced  by  a  white  (in  time  and  space)  additive  Gaus¬ 
sian  noise.  In  our  case,  however,  no  peak  appears  in  the  spec¬ 
tral  density  (a2  )  at  a  given  wavenumber  m  for  an  intermediate 
level  of  noise. 

Our  results  show  the  coherent  behavior  of  a  group  of 
modes  when  amplitudes  of  many  modes  in  the  group  ex¬ 
ceed  a  threshold.  Low-  and  high-order  modes  are  separately 
grouped.  For  quite  large  amplitudes  our  results  demonstrate 
the  driver-response  relationship  for  which  phase  of  modes 
are  locked. 

To  examine  this  mechanism  we  have  linearized  govern¬ 
ing  Eqs.  (6-8),  and  calculate  the  growth  of  perturbations 
excited  by  stochastic  forcing  Eq.  (10)  in  this  case.  The  re¬ 
sults  of  these  calculations  are  summarized  as  follows.  The 
power  growth  of  perturbations  in  the  transient  regime  is  ob¬ 
served  up  to  day-20  and,  therefore,  cannot  be  caused  by  the 
phase  locking  mechanism.  “Super-exponential”  growth  of 
perturbations  after  day-25  disappears  when  nonlinear  (iner¬ 
tial)  terms  are  removed  from  the  governing  equations. 

The  power  growth  of  perturbations  in  the  transient  regime 
can  be  explained  using  results  obtained  by  Seki  et  al.  (1993). 
They  pointed  out  that  for  a  linear  dynamical  system  forced 
by  a  Gaussian  white  noise  the  mean  deviation  of  perturbation 
amplitudes  can  grow  along  a  power  law  up  to  a  time  scale  7); 
defined  by  inverse  of  friction  coefficient.  For  large  values 
of  friction  coefficients  such  a  power-law  behavior  disappears 
because  7),— >0. 

The  effective  dissipation  in  model  (6-8)  depends  on  the 
structure  of  the  reference  flow  and  shape  of  the  spectrum  of 
Rm .  When  forcing  becomes  stronger  or  is  inhomogeneous, 
the  center  mass  of  energy  spectrum  shifts  to  high  wavenum¬ 
ber  domain  and,  in  general,  the  dissipation  of  perturbations 
reduces  because  the  nonlinear  bottom  friction  is  most  effec¬ 
tive  for  largest-scale  perturbations  as  it  has  been  checked  nu¬ 
merically.  That  results  into  appearance  of  power-law  behav¬ 
ior  for  variance  of  perturbations.  In  this  case  power  exponent 
depends  on  the  structure  of  time-dependent  reference  flow, 
and  unfortunately,  cannot  be  analytically  calculated.  Seki 
et  al.  (1993)  calculated  the  power  exponents  for  two  simple 
stochastic  dynamical  systems  only.  However  explicit  corre¬ 
lation  between  level  of  model  dissipation  and  existence  of  a 
power-law  behavior  of  perturbations  are  clearly  observed  in 
our  numerical  experiments. 

After  25  day  integration  noiseless  model  (6-8)  reaches  a 
spin  up  when  a  solution  oscillated  with  a  period  of  120  days 
for  the  dominant  low-order  modes  within  any  1000-day  time 
interval.  Oscillating  modes  weakly  interact  one  with  another 
due  to  nonlinear  (inertial  and  frictional)  terms  in  the  gov¬ 
erning  equations.  Effects  of  external  noise  on  these  modes 
lead  to  coherent  behavior  of  modes  that  seems  to  be  simi- 
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lar  to  the  phase  locking  of  the  modes,  which  is  understood 
in  a  statistical  sense,  as  the  existence  of  a  preferred  value  of 
the  phase  difference  between  individual  stochastically  forced 
modes  with  weak  interactions  among  them. 

Non-steady  nature  of  the  reference  flow  is  essential  for  ex¬ 
istence  of  super-exponential  regime  of  perturbation  growth. 
This  was  checked  by  numerical  modeling.  The  oscillations 
are  smoothed  and  the  reference  flow  becomes  steady  when 
the  drag  coefficient  is  2-3  times  as  much.  Neither  “super  ex¬ 
ponential’’  growth  of  perturbations  nor  coherent  behavior  of 
modes  is  observed  for  such  levels  of  model  dissipation. 

Note  that  high-order  modes  of  the  reference  flow  have  a 
dominant  oscillation  period  of  about  30  days,  not  120  days 
as  low-order  modes.  Since  the  first  twenty  modes  contain  up 
to  90%  of  the  kinetic  energy  of  circulation,  120-day  oscilla¬ 
tions  dominate  the  flow  and  mask  faster  motions.  However, 
existence  of  the  second  dominant  period  leads  to  coherence 
in  behavior  of  high-order  modes,  which  is  different  than  that 
low-order  modes  demonstrate. 

The  coherence  in  behavior  of  modes  can  be  explained 
by  a  number  of  mechanisms,  such  as  synchronization  of 
weakly  coupled  oscillators  (Pikovsky  et  al.,  2000),  modu¬ 
lation  (Landa,  1996)  and  others.  In  practice  these  mech¬ 
anisms  cannot  be  selected  using  only  a  driver-response  re¬ 
lationship  or  the  cross-spectral  technique  without  a  simple 
physical  model.  Unfortunately,  this  is  a  great  problem  to 
develop  a  model,  which  would  adequately  describe  nonlin¬ 
ear  oceanic  flow  dynamics  with  many  degrees  of  freedom. 
Therefore  herein  we  are  not  able  to  select  one  of  the  physical 
mechanisms  discussed  in  modern  literature  for  explanation 
of  the  coherent  behavior  of  modes  observed  in  our  numerical 
experiments  for  finite-amplitude  perturbations.  However,  un¬ 
doubtedly  the  observed  coherent  behavior  is  due  to  nonlinear 
interactions  among  modes  and  oscillations  of  the  reference 
flow. 

7  Finite-amplitude-induced  transition  in  predictability 
skill 

Our  computations  have  shown  that  for  strong  winds  model 
predictability  demonstrates  stronger  sensitivity  to  amplitudes 
of  perturbations  induced  by  the  stochastic  wind  than  to  the 
choice  of  spatial  structure  function  in  Eq.  (10).  Intuitively, 
larger-amplitude  perturbations  should  cause  faster  decay  of 
model  predictability,  what  was  confirmed  by  our  numerical 
experiments  for  the  mean  predictability  time,  r-mean  mono- 
tonically  reduced  as  er2  increases  (not  shown). 

Furthermore,  our  experiments  have  also  shown  that  the 
collective  behavior  of  finite-amplitude  perturbations  may 
cause  sudden  changes  (bifurcations)  in  the  high-order  statis¬ 
tics  of  predictability  time.  This  effect  is  clearly  observed 
when  perturbation  amplitudes  exceeded  some  threshold,  af¬ 
ter  which  the  global  correlations  among  the  perturbations 
with  different  scales  dominated  the  PE  characteristics. 


Fig.  12.  r-statistics  for  different  values  of  e~  and  a2,  (a)  r-mean, 
(b)  r-variance,  (c)  r-skewness  and  (d)  r-kurtosis.  Triangles,  aster¬ 
isks  and  circles  correspond  to  cr2=2.0,  1.0,  and  0.5. 

We  have  called  such  bifurcations  as  the  “finite-amplitude 
induced  phase  (non  thermodynamic)  transitions  in  model 
predictability”.  They  are  detected  using  statistics  of  IPT, 
such  as  r-variance,  r-skewness  and  r-kurtosis,  but  other 
measures  of  model  predictability  are  also  applicable. 

Finite-amplitude  phase  transitions  should  easily  be  de¬ 
tected  for  any  hydrodynamic  model  using  the  dependence  of 
IPT  statistics  on  e2  because  the  value  of  the  tolerance  level 
limits  the  maximum  amplitude  of  perturbations  existing  in 
the  model.  As  an  example,  detection  of  the  phase  transition 
for  model  (6-8)  is  described  below. 

For  small  tolerance  (e2<5.0x  10-3)  model  predictability 
is  low:  r-mean  does  not  exceed  20  days  (Fig.  12a)  and  r- 
variance  is  quite  large,  up  to  25-27  days2  (Fig.  12b).  Ad¬ 
ditionally,  a  large  negative  skewness  about  —0.5  (Fig.  12c) 
indicates  that  the  IPT  distribution  has  a  tail  stretching  into 
domain  of  small  prediction  times. 

The  mean  IPT  monotonically  grows  as  e2  increases 
(Fig.  12a).  Although  here,  the  IPT  grows  with  various  rates 
for  different  tolerance  levels,  no  bifurcations  are  observed  in 
this  figure.  In  contrast  to  r-mean,  the  value  of  r-variance 
suddenly  changes  when  e2  becomes  larger  than  5.0xl0-3 
(this  value  is  taken  as  a  threshold).  The  variance,  which  was 
quite  large  (about  25-27  day2)  for  small  tolerances,  suddenly 
reduces  to  5  day2  when  e2  crosses  the  threshold  (Fig.  12b). 

Both  r-skewness  and  r-kurtosis  also  change  considerably 
(Figs.  12c,  d).  They  converge  asymptotically  to  1.0  and  4.7, 
respectively,  as  e2  increases.  Positive  skewness  corresponds 
to  asymmetric  r-PDF  shapes  with  a  tail  stretching  into  large 
prediction  times.  The  large  kurtosis  (much  larger  than  3)  in¬ 
dicate  that  PDF  is  highly  non-Gaussian. 

Typical  r-PDFs  computed  before  and  after  the  phase  tran¬ 
sition,  are  given  in  Figs.  13a  and  b,  respectively.  Comparing 
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x  (days)  x  (days) 


Fig.  13.  Histograms  of  IPT  computed  for  <r2=1.0. 
(a)  e2=1.0x  10-2,  (b)  e2=0.1. 


them  one  to  another,  we  find  that  the  phase  transition  causes 
the  PDF  tail  stretching  into  domain  of  large  prediction  times. 
This  tail  is  formed  by  rare  predictions  of  duration  up  to  50 
days.  Smaller  variance  and  larger  positive  skewness  show 
that  model  predictability  was  considerably  enhanced  after  the 
phase  transition. 

Two  physical  mechanisms  are  responsible  for  the  phase 
transition.  First,  it  is  clearly  from  Fig.  12a  that  the  phase 
transition  exists  when  the  mean  r-IPT  is  not  less  than  30 
days.  During  this  time  period  the  model  reaches  the  spin  up 
state.  Different  stability  properties  in  the  phase  space  near 
and  far  away  from  the  point  are  caused  by  the  inhomogene¬ 
ity  of  model  phase  space.  That  leads  to  different  statistics  of 
PE  before  and  after  the  phase  transition.  Second,  even  when 
perturbation  amplitudes  are  several  percents  of  the  reference 
solution,  a  coherent  behavior  of  different  scales  are  clearly 
observed.  In  this  case  r-variance  reduces  due  to  strong  cor¬ 
relations  among  the  scales  (Kravtsov,  1993). 


8  Weibullian  statistics  of  IPT 

Our  computations  have  shown  that  stochastic  forcing 
Eq.  (10),  in  general,  induces  highly  non-Gaussian  r-PDFs  for 
finite-amplitude  PEs.  The  following  question  arises:  what 
kind  of  statistics  can  be  used  to  represent  such  r-PDFs?  If 
appropriate  distribution  function  is  found,  it  would  be  pos¬ 
sible  to  identify  the  ensemble  generated  PDFs  from  limited 
observation  series  and  small  forecast  ensembles,  and  in  turn 
to  estimate  the  model  predictability  horizon  (i.e.  maximum 
predictability  time  reached  for  the  given  model  and  wind  un¬ 
certainty  (Kravtsov,  1993)). 

Nonlin.  Processes  Geophys.,  14,  1-16,  2007 


We  apply  the  three-parameter  Weibull  statistics  with  dis¬ 
tribution  f{r)  and  cumulative  distribution  Pit) 


for  the  analysis  of  r.  Here,  ??,  y,  and  ft  are  scale,  shape, 
and  location  parameters  (von  Storch  and  Zwiers,  1999).  The 
following  original  algorithm  is  developed  to  estimate  the  pa¬ 
rameters  of  distribution  (27)  from  an  IPT  ensemble  sampling. 

Closeness  of  two  inverses  of  r-CDFs:  X  (P)(“real”  value) 
and  Xq(P)  (the  first  guess),  may  be  estimated  by  the 
Kullback-Leibler  distance  Q  (White,  1994): 

l 

Q  =  j  X(P)\n[X(P)/X0(P)]clP.  (29a) 

o 

Then,  X(P)  is  the  solution  of  the  following  variation  prob¬ 
lem  (Kapur  and  Kesavan,  1992) 

— »  min .  (29b) 

The  Kullback-Leibler  distance  Q  should  be  a  subject  to  ad¬ 
ditional  constraints  from  the  following  condition:  the  proba¬ 
bility  weighted  moments  computed  from  the  ensemble  sam¬ 
pling  (ai,  a 2,  and  aj)  and  theoretically  (dq,  cG.and y)  from 
Eq.  (4)  must  coincide.  This  condition  is  accounted  for  in 
Eq.  (29)  through  additional  constrain  as 

I 

n  =  y  X(P)ln[X(P)/X0(P)]dP  +  xi  («t-at) 
o 

+X2  (a 2  -  a2)  +  X3  (“3  -  a 3)  min  .  (30) 

where  xi,  X2-  and  X3  are  Lagrange  multipliers. 

Functional  Eq.  (30)  is  minimized  with  respect  to  X(P). 
The  solution  of  minimization  problem  (Eq.  30)  is  written  as 

X(P)  =  Xq(P)  exp  (-Xi P  ~  X2 P2  ~  X3^3)  -  (31) 

For  details  see  Kapur  and  Kesavan  (1992).  Then,  the  La¬ 
grange  multipliers  /  \ ,  X2.  and  X3  are  determined  by  the 
quasi-Newton  iteration  method  as  a  solution  of  nonlinear 
least-square  problem  resulting  from  Eqs.  (30)  and  (31). 

Our  computations  show  that  (a)  the  method  discussed 
above  is  robust  relative  to  sampling  error  if  only  few  mo¬ 
ments  are  used  as  constrains,  and  (b)  Lagrange  multipliers 
are  estimated  within  10-12  iterations  only.  In  general  case, 
when  sampling  error  is  considerable  and  more  moments  are 
required  in  Eq.  (30),  the  non-linear  least-squares  minimiza¬ 
tion  problem  is  solved  through  the  Levenberg-Marquardt  it¬ 
erative  method  (Engl  et  ah,  1996). 
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We  did  not  find  difference  between  the  distribution  func¬ 
tion  calculated  by  a  non-parametrical  technique  based  on  the 
Epanichenikov’s  kernel  and  the  bootstrap  re-sampling  pro¬ 
cedure  (Good,  2001)  directly  from  ensemble  sampling,  and 
appropriate  Weibull  counterpart  obtained  by  the  method  (29- 
31),  at  least  at  95%  confidence  level. 

This  conclusion  is  illustrated  for  the  ensemble  sampling 
shown  in  Fig.  14a.  The  r-CDFs  computed  by  the  non- 
parametrical  technique  (solid  curve)  and  our  method  (black 
dots)  are  compared  in  Fig.  14b.  Differences  between  them 
are  negligible.  The  parameters  p,  y,  and  /I  are  equal  to 
37.1xl0_1  days,  30.0  days  and  16. 7 x  10“',  respectively. 

Parameter  /)  affects  the  length  of  the  PDF  tail  formed  by 
rare  forecasts,  which  are  longer  than  the  mean  ensemble  fore¬ 
cast  (r).  Small  f3  indicates  enhanced  probability  for  real¬ 
ization  of  abnormal  long  (in  our  case  up  to  50  days)  model 
forecasts. 

9  Model  predictability  horizon 

Asymptotic  behavior  of  r-CDF  as  r-*o o  determines  the 
predictability  horizon  of  the  model,  i.e.  the  maximum  pre¬ 
dictability  time  of  an  individual  forecasting  for  the  given 
model  and  statistics  of  wind  perturbations  (Kravtsov,  1993). 

Accordingly  to  Eq.  (28)  the  model  predictability  horizon 
is  calculated  by 

Thor  =  y  +  i]\-\nP*]l/p,  (32) 

where  P*  is  the  probability  that  thor  will  be  achieved  in  an 
individual  forecasting.  For  fixed  P*  Eq.  (32)  shows  a  slow 
power-law  growth  of  the  predictability  horizon  with  the  de¬ 
crease  of  the  shape  parameter. 

Fet  us  estimate  the  predictability  horizon  for  the  example 
discussed  above.  Substitution  of  the  distribution  parameters 
obtained  above  into  Eq.  (32)  leads  to 

^hor  ^  40.2  days,  45 .4  days  and  50.5  days  (33) 

for  P*=1.0xl0“2,  1.0xl0“3  and  l.OxlO-4,  respectively. 
These  estimations  demonstrate  that  for  the  chosen  values  of 
s 2  and  a,  the  model  predictability  horizon  is  limited  to  50 
days,  and  any  individual  forecasting,  which  is  longer  than  50 
days,  is  unlikely. 

10  Conclusions 

A  simple  shallow-water  model  was  used  to  understand  sen¬ 
sitivity  and  predictability  of  ocean  models  with  inaccurate 
wind  forcing.  This  model  used  a  highly  idealized  represen¬ 
tation  of  ocean  dynamics  and  did  not  simulate  the  redistri¬ 
bution  of  PE  between  barotropic  and  baroclinic  dynamics  as 
high-resolution  ocean  models.  However,  due  to  the  small 
number  of  degrees  of  freedom  of  the  model  (only  462  vari¬ 
ables),  distribution  functions  for  predictability  scale  and  its 


Fig.  14.  Identification  of  r-CDF.  (a)  IPT  histogram  for  a  103  term 
ensemble;  (b)  CDF  computed  directly  from  the  ensemble  (solid 
curve)  and  using  the  developed  method  (black  dots),  y=30  days. 

high-order  moments  were  computed  for  a  large  number  of 
ensemble  realizations  (up  to  50000).  This  guaranteed  re¬ 
duced  sampling  error  and  robustness  in  estimating  the  PE 
statistics. 

Similar  analysis  is  difficult  to  undertake  in  full-scale  nu¬ 
merical  forecast  ocean  models  due  to  limited  computer  re¬ 
sources.  Generally,  the  full-scale  models  produce  small  en¬ 
semble  samples  and  therefore,  cannot  resolve  the  full  com¬ 
plexity  of  the  PE.  The  idealized  model  revealed  trends  in  PE 
behavior  and  the  model  reconstructed  PE  statistics  with  min¬ 
imum  distortion.  These  statistics  can  be  used  for  the  analysis 
of  the  smaller  ensemble  samples  from  the  full-scale  models. 

Baroclinic  high-resolution  ocean  models  should  be  used  to 
examine  the  following  trends  obtained  in  the  present  study. 

The  predictability  time  for  small  perturbations  may  be 
much  larger  than  the  inverse  of  the  leading  local  character¬ 
istic  exponent.  The  shallow  water  model  showed  that  it  is 
possible  to  have  non-trivial  time  evolutions  of  small  (but  fi¬ 
nite)  perturbations  and  that  their  growth  could  be  fitted  by 
power  laws  although  the  perturbations  were  actually  ampli¬ 
fied  by  the  background  flow.  The  power  growth  for  all  or  a 
portion  of  scales  was  determined  by  the  cumulative  effects  of 
multiple  characteristic  times  . 

The  expected  growth  of  error  and  decay  of  skill  occurs 
most  rapidly  for  smaller  scales  and,  with  time,  expands  to 
larger  scales.  One  of  the  main  features  of  the  ocean  is  the  ex¬ 
istence  of  strongly  interacting  spatial  scales,  which  raises  the 
possibility  of  different  behavior  of  the  PE  at  different  scales 
of  motion.  It  is  traditionally  assumed  that  small  scales  are 
less  predictable  than  larger  scales.  This  picture  was  drawn 
by  Lorenz  (1969)  for  perfect  model  scenario,  and  then  trans¬ 
ferred  to  the  analysis  of  forecast  error  as  a  function  of  spatial 
scale  in  operational  atmospheric  (see,  for  example  Dale  her 
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and  Kalnay,  1987)  and  oceanographic  models  (Brasseur  et 
al.,  1996,  among  others). 

Our  simplified  model  predicts  existence  of  two  additional 
predictability  regimes  for  imperfect  models.  For  large  scale 
stochastic  winds,  the  perturbations  rapidly  grew  to  horizontal 
scales  comparable  to  those  of  the  reference  state.  In  contrast, 
smaller  scale  perturbations  excited  by  the  wind  were  subject 
to  strong  viscous  damping.  Therefore,  the  predictability  of 
wind-driven  circulation  was  less  affected  by  model  uncer¬ 
tainties  acting  at  small  scales  than  at  larger  scales.  Mahade- 
van  et  al.  (2001),  using  a  quasi-geostrophic  ocean  circula¬ 
tion  model  for  perfect-model  twin  experiments,  found  that 
such  a  scenario  was  favourable  for  weakly  aperiodic,  peri¬ 
odic,  and  stationary  circulation  regimes  when  the  mesoscale 
energy  content  was  relatively  low. 

For  stochastic  wind  that  is  limited  to  scales  smaller  than 
those  occupied  by  large-scale  flow,  perturbations  on  differ¬ 
ent  scales  may  grow  coherently  due  to  interactions  among 
them.  The  coherent  growth  of  perturbations  has  been  iden¬ 
tified  on  different  scales  at  various  stages  of  PE  evolution. 
Coherent  behavior  of  PE  can  be  found  in  full-scale  oper¬ 
ational  atmospheric  models  (Boer,  2003,  as  an  example), 
and  quasi-geostrophic  models  (Vannitsem  and  Nicolis,  1997; 
McWilliams  and  Chow,  1981,  among  others).  McWilliams 
and  Chow  (1981)  demonstrated  that  for  a  simple  three-level 
quasi-geostrophic  model,  all  scales  of  motion  exhibited  a 
similar  growth  rate  after  a  short  transient  phase.  The  coherent 
growth  of  PE  on  different  scales  can  be  found  in  imperfect 
quasi-geostrophic  models  too  (for  example,  see  Vannitsem, 
2006).  These  and  other  results  (not  discussed  here)  indicate 
that  the  coherent  growth  of  PE  may  be  caused  by  model  in¬ 
dependent  (universal)  mechanisms,  which  require  further  in¬ 
vestigation. 

The  present  study  introduced  a  new  statistics  (Weibull) 
for  finite-amplitude  PE  and  suggested  a  practical  way  for  its 
identification  through  the  probability  weighted  moments  and 
a  variation  principle.  This  extremum  statistics  is  often  ob¬ 
served  to  arise  in  finite  sized,  multi-body  systems,  exhibiting 
correlation  over  a  broad  range  of  scales,  leading  to  emergent 
phenomenology,  such  as  self-similarity  and  in  some  cases 
fractional  dimensions  (Boffetta  et  al.,  2002).  A  universal  ap¬ 
proach  to  extract  extremum  statistics  from  short-  and  inter¬ 
mediate  marine  forecasts  was  suggested.  Possible  general¬ 
ization  of  the  approach  for  small  forecast  ensembles  will  be 
discussed  in  a  separate  paper. 


Appendix  A 


where  Zm  and  Qm  are  linear  and  nonlinear  (inertial  and  fric¬ 
tional)  terms,  respectively  (for  details  see  Pedlosky,  1987), 
and 


r m  —  Ym  Z  (bm  sin  0  Cm  cos  0 )  . 


(A2) 


Ym  —  —CdH  Va2,  bm=  II  i — ifmdx\dx2, 
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The  stochastic  processes  z(t )  and  9(t)  are  defined  in  Sect.  4. 
The  variance  of  PE  <a2  >  satisfies  the  obvious  equation 


d<ar„> 


dt 


— 2  (  <  Zm  Clm  >  -f-  <  Q  m  Clm  >  ~t“  <  I'm  Clm  >  ) , 


(A4) 


with  the  wind  error  source  term 


Rm  =  Ym  [bm  (z  sin  9am)  +  cm  (z cos  9am >] .  (A5) 


Since  /(z,  9)—f(z)f(9),  average  of  Eq.  (A5)  can  be  divided 
into  two  steps.  (1)  rmam  is  averaged  over  z.  (2)  The  obtained 
function  <rmam>z  is  averaged  over  9. 

The  correlation  function  <zam>z  is  analytically  calcu¬ 
lated  using  the  cumulant  decomposition  (Klyatskin,  2005): 


{z(t)am(t))Wr 


dti...dtsKs+i(t,  t\ ,  ...,  ts) 


I  Ssam(t)  \ 
\Sz(h)...8z(ts)lz ‘ 


(A6) 


Here,  Ks(t\, ...,  ts)  is  the  s-th  order  cumulant  of  the  noise  z. 
The  notation  jZ  denoted  the  functional  derivative. 

For  a  stationary  8— correlated  noise  z  one  obtains 


Ksih,  ...,  ts )  =  ks(t\)8(ti  -  t2)...S(ts  -  ts- 1),  (A7) 


where  ks(t i)  are  the  intensity  coefficients  (Stratonovich, 
1963). 

Substituting  Eq.  (A7)  into  Eq.  (A6)  yields 
(z(t)8am(t))z  =  (A8) 


From  Eq.  (Al)  we  obtain  for  1 


/  (t)\ 

\  Sz(t)  lz 


Ym  (bm  cos  9  +  cm  sin  9) , 


(A9) 


Wind  error  source  term 


and  for  s> 2 


Using  the  quasi-geostrophic  approximation,  governing 
Eqs.  (6-8)  can  be  re-written  in  the  spectral  form 

d  ci  fjj 

.  —  Zmicik)  A  Qm  (t7/,  Uk)  T  rm  (t),  (Al) 


I8sam(t)\ 

\  Sz(t)s  lz 


(A10) 


Taking  into  consideration  that  for  exponential  distribution 
function  the  s-order  cumulant  is  calculated  as  ks—2-(s  —  1)! 
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(Zelen  and  Severo,  1972),  and  averaging  Eq.  (A9)  over  the 
stochastic  process  6  ( t )  we  find 

=  vlibi  +  4),  (All) 

The  coefficients  bm  and  cm  depend  only  on  the  spatial  in¬ 
homogeneity  of  stochastic  wind  forcing.  Spectrum  of  wind 
error  term  Rm  for  different  structure  functions  G  is  shown  in 
Figs.  9a,  b,  c. 
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